Characterization of Coastal Systems#

Welcome to the first notebook exercise of Coastal Systems (TU Delft, MSc Coastal Engineering)! This is the first year that we will experiment with notebooks in this course. With these notebooks we hope to provide you with interactive course material that helps you better understand the processes and concepts that we teach in this course. Please let us know how you find the notebooks - we appreciate your feedback!

Chapter 2 of Coastal Dynamics Open Textbook describes the large geographical variation of coasts across the world. It explains how the coasts that we have today are shaped by both present-day processes and processes millions years ago. It distinguishes between three different order of features, which are associated to different orders of of time. In this notebook we will look at coastal systems at these different orders of scale.

Import libraries that we use for our analysis#

In the two cells below we import the libraries that we need for the analysis. We also set some path settings to load the data and source code. For example, in the cell below we add the src directory to the system path, which allows us to import generic functions from ../../src/coastpy.

import os
import pathlib
import sys

# make coastpy library importable by appending root to path
cwd = pathlib.Path().resolve()
proj_dir = cwd.parent.parent  # this is the root of the CoastalCodeBook
sys.path.append(str(proj_dir / "src"))

import warnings

os.environ["USE_PYGEOS"] = "0"  # use shapely 2.0 instead of PYGEOS and silence warnings

import colorcet as cc
import dask.dataframe as dd
import geopandas as gpd
import holoviews as hv
import hvplot.pandas  # noqa: API import
import hvplot.xarray  # noqa: API import
import ipyleaflet
import numpy as np
import pandas as pd
import panel as pn
from geoviews import tile_sources as gvts
from ipyleaflet import Map, Marker, ScaleControl, basemaps
from ipywidgets import HTML

pn.extension()

from coastpy.geometries import geo_bbox

# definitly not a best practice, but a workaround to avoid many warnings triggered by shapely 2.0 release
warnings.filterwarnings("ignore", category=RuntimeWarning)
warnings.filterwarnings("ignore", category=UserWarning)

DATA_DIR = proj_dir / "data"
coastal_systems_fp = DATA_DIR / "01_coastal_systems.gpkg"

Exercise 1: Plate tectonics#

In this notebook we will start with the broadest (or first order) features of the coast that cover large geographical distances (thousands of kilometres) and are linked to the long-term geological process of plate tectonics. We will do so by using earthquake data from the USGS. The dataset we load contains a sample (10%) of observed eartquakes between Jan 2000 and Dec 2018. Why earthquake data? Earthquake data reveals geologists the mysteries of the deep, but also for coastal researchers the data is insightful. Let’s first load the data by running the next cells.

Load the earthquake data#

We load the data (tabular data including geometries) and index the columns to only keep the data in memory that we actually need. In total the dataset contains 2.2 million earthquakes, but here we use a sample (10%), so the data contains approx. 220k eartquake entries. If you find that the interactive panel responds slow to how you adjust the widgets, please consider to take another sample. You can do so by uncommenting the sample line in the next cell. So if you set frac=0.1 you have a dataframe with approx. 22k unique earthquakes over the world.

WEB_MERCATOR_LIMITS = (
    -20037508.342789244,
    20037508.342789244,
)  # max polar latitudes that can be handled in World Mercator

df = (
    dd.read_parquet(DATA_DIR / "01_earthquakes_sample.parquet")
    .sample(
        frac=0.1
    )  # uncomment this line if loading the data takes too long on your computer
    .set_index("time")
    .compute()
    .tz_localize(None)
    .sort_index()
)


# To save memory we drop most of the columns. Also we drop the polar latitudes that cannot be displayed in the web mercator projection.
df = df[["mag", "depth", "latitude", "longitude", "place", "type"]][
    df["northing"] < WEB_MERCATOR_LIMITS[1]
]
df.head()

	mag	depth	latitude	longitude	place	type
time
2000-01-02 18:48:55.600	3.90	32.000	35.480	22.510	central Mediterranean Sea	earthquake
2000-01-03 03:43:08.900	2.05	7.399	34.725	-116.301	13km W of Ludlow, California	earthquake
2000-01-03 05:47:06.280	1.40	1.602	34.272	-116.446	18km N of Yucca Valley, California	earthquake
2000-01-04 08:13:20.000	2.00	11.000	38.100	-90.480	7km ESE of De Soto, Missouri	earthquake
2000-01-04 20:48:47.450	1.70	9.699	34.670	-116.292	13km WSW of Ludlow, California	earthquake

Visualization of the earthquake data#

To explore the data we use visualization tools from the Holoviz project that makes high-level tools to simplify visualization in Python. In the next cell we enable the interactive mode on the data dataframe, create widgets to explore the data and filter the dataframe accordingly. To explore the eartquake data we create an overlay of the eartquakes on a tileset of ESRI Imagery. Please note that the code in the next cell will only do the computations and store the result in an object called panel. To actually see the results you have to run one more cell; the one that calls this object panel.

pn.extension()
title_bar = pn.Row(
    pn.pane.Markdown(
        "## Exercise 1: Plate tectonics & first-order coastal features",
        style={"color": "black"},
        width=800,
        sizing_mode="fixed",
        margin=(10, 5, 10, 15),
    ),
    pn.Spacer(),
)

# define widgets that can be used to index the data
magnitude_slider = pn.widgets.RangeSlider(
    name="Earthquake magnitude [Richter]", start=0.1, end=10
)
depth_slider = pn.widgets.RangeSlider(name="Earthquake depth [km]", start=0.1, end=650)
date_slider = pn.widgets.DateRangeSlider(
    name="Date", start=df.index[0], end=df.index[-1]
)
column_types = pn.widgets.Select(options=["mag", "depth"])


@pn.depends(
    magnitude_slider.param.value_start,
    magnitude_slider.param.value_end,
    depth_slider.param.value_start,
    depth_slider.param.value_end,
    date_slider.param.value_start,
    date_slider.param.value_end,
    column_types.param.value,
)
def plot_earthquake_panel(
    magnitude_start,
    magnitude_end,
    depth_start,
    depth_end,
    date_start,
    date_end,
    column_type,
):
    panel = df[(df.mag > magnitude_start) & (df.mag < magnitude_end)]

    panel = df[
        (df.mag > magnitude_start)
        & (df.mag < magnitude_end)
        & (df.depth > depth_start)
        & (df.depth < depth_end)
        & (df.index >= date_start)
        & (df.index <= date_end)
    ]
    # inverted fire colormap from colorcet
    cmap = cc.CET_L4[::-1]
    colorbar_labels = {"mag": "Magnitude [Richter]", "depth": "Earthquake depth [km]"}

    return panel.hvplot.points(
        x="longitude",
        y="latitude",
        geo=True,
        color=column_type,
        global_extent=True,
        tiles="ESRI",
        frame_width=900,
        ylabel="Latitude [deg]",
        xlabel="Longitude [deg]",
        cmap=cmap,
        tools=["tap"],
        hover_cols=["place", "time"],
        logz=True,
        clim=(1, None),
        clabel=colorbar_labels[column_type],
    )


cmap = cc.CET_L4[::-1]  # inverted fire colormap from colorcet

earthquake_panel = pn.Column(
    title_bar,
    pn.Row(column_types),
    pn.Row(pn.Column(magnitude_slider, depth_slider)),
    pn.Row(depth_slider),
    pn.Row(date_slider),
    pn.Row(plot_earthquake_panel),
)

If the visualization is too slow, please follow the instructions in loading the data for taking a sample.

After running the cell below you will have a panel with several widgets to index the eartquake data; by magnitude, depth and time, while the colors on the map show either the magintude or the depth of the earthquakes.

earthquake_panel

Explore the earthquake data & questions#

How do the earthquake magnitude and earthquake depth relate to the coasts that we see? (Hint: See Figure 2.3 in the textbook and consider how deep under the ground the plates are moving. Extra hint: How do earthquake magnitude and depth differ for convergent and divergent plate boundaries?)
Earthquake data support one of the most fundamental processes in the geology: plate tectonics. Although plate tectonics is a relatively slow process that acts on the geological time scale, it has had an enormous impact on the formation of coastlines and determines the broadest features of the coast. What are two important inherited aspects of this process? (Hint: see Figure 2.10 and Sec. 2.3.3 in the textbook.)
In 1971 Inman, D. L. & Nordstrom, C. E. used plate tectonics to classify the coast. Explain the classification that they introduced. What are the three different classes that they distinguish? How do they match with the earthquake data as you can explore in the panel?
Can you identify or predict areas around the world where you will find the coasts that are distinguished by Inman, D. L. & Nordstrom, C. E.? For instance, what kind of coasts do you have in Chili? And how are they different to the east coast of the USA? And what is characteristic about the East China sea?
Inman, D. L. & Nordstrom (1971) further distinguish Afro-trailing-edge coasts and Amero-trailing-edge coasts based on differences in sediment supplies. What is the main cause of these differences in sediment supply? And how do you expect the differences in sediment input to show in the coastal geomorphology?

Exercise 2: Process-based coastal classification#

In the section part of this notebook we will explore several coastal systems around the world considering the second and third order scale of features. In chapter two of the Coastal Dynamics open textbook it is explained how coastal systems can be classified according to the processes that characterize these systems. For example, one of the figures (below) shows how the relative influence of fluvial, wave and tidal processes influences the shape of coastal features. The idea of this exercise is that you identify the signatures of the processes that are introduced in chapter 2 in several coastal systems around the world.

The coastal systems data#

In the cell below we define a small plot function that generates a ESRI World Imagery basemap given a longitude, latitude, zoom level and name. Also, a small sheet of coastal systems around the world is loaded into geopandas, a Python library for geospatial tabular data. In the cells afterwards we sample this dataframe and show the coastal system on a map. Since the sample is random you might encounter the same coastal system multiple times; then you can just run the cell again to get another ‘coastal draw’.

coastal_systems = gpd.read_file(coastal_systems_fp)

pn.extension()
title_bar = pn.Row(
    pn.pane.Markdown(
        "## Exercise 2: Coastal system characterization",
        style={"color": "black"},
        width=800,
        sizing_mode="fixed",
        margin=(10, 5, 10, 15),
    ),
    pn.Spacer(),
)

options = coastal_systems.name.to_list()
coastal_systems_slider = pn.widgets.Select(
    name="Coastal system", options=options, value=np.random.choice(options)
)


@pn.depends(coastal_systems_slider.param.value)
def plot_coastal_system(name):
    system = coastal_systems.loc[coastal_systems["name"] == name].copy()
    west, south, east, north = system[
        ["west", "south", "east", "north"]
    ].values.flatten()

    return system.hvplot.points(
        x="lon",
        y="lat",
        geo=True,
        color="red",
        alpha=0,
        xlim=(west, east),
        ylim=(south, north),
        tiles="ESRI",
        frame_width=1100,
        ylabel="Latitude [deg]",
        xlabel="Longitude [deg]",
    )


app = pn.Column(
    title_bar,
    pn.Row(coastal_systems_slider),
    pn.Row(plot_coastal_system),
)
app

Explore the coastal systems#

While sampling over a range of coastal systems, try to answer the following questions.

Find and compare a heavily engineered river-dominated delta and a more natural river-dominated delta
Compare the scale of the biggest and smallest tidal basin in the dataset
Find the estuarine and deltaic systems with a spit
Compare and contrast wave-dominated deltas with high and low sediment supply. How can you tell?
Find a tidal estuary with large fine (muddy) sediment supply, then find one with a large coarse (sandy) sediment supply. How can you tell the difference?
The eastern and western tips of the Dutch and German Wadden Islands are very different beach ridge environments. How might differences in sediment supply explain this? Where is the sediment coming from?
The Dune du Pilat in France is one of the world’s largest coastal sand dunes (it is also one of the coolest and you should definitely visit if you get the chance!). Why is it located on the east side of Arcachon Inlet and not the west?
Look at the northern Jiangsu coast in China. What might explain the limited sediment supply in this location?
Find an estuary or tidal bay with extensive intertidal flats. Do you see salt marshes or mangrove forests nearby? Why or why not?
Find an inlet with jetties. How might this affect the way it evolves?
Find a delta/estuary/inlet whose shape is constrained by the presence of rocky coastal features.
The Albufeira Lagoon in Portugal opens and closes seasonally. In the image shown, is it open or closed? When and how might it open or close?
Find examples of heavily urbanized estuaries. How might these human interventions influence the natural processes there?
Based on these satellite images, which is the most beautiful site? Taking a moment to appreciate the beauty of these natural systems is an important part of your job as coastal engineers.

from coastpy.maps import plot_esri_basemap

m = plot_esri_basemap(-9.181333, 38.510379, 13, "Lagoa de Albufeira")
m

bbox = geo_bbox(m.west, m.south, m.east, m.north)

bbox.explore()

Make this Notebook Trusted to load map: File -> Trust Notebook

m.west, m.south, m.east, m.north

(0.0, 0.0, 0.0, 0.0)

import shapely
albufeira = {
    "name": "Lagoa de Albofeira",
    "lon": -9.181333,
    "lat": 38.510379,
    "zoom": 13,
    "west": -9.276924133300783,
    "south": 38.483560395392516,
    "east": -9.085865020751955,
    "north": 38.53729004998249,
    "geometry": shapely.geometry.point.Point([-9.181333, 38.510379])
}
albufeira = gpd.GeoDataFrame([albufeira], crs=4326)

coastal_systems = pd.concat([albufeira, coastal_systems])
coastal_systems = coastal_systems[coastal_systems["name"]!="Albufeira"]

coastal_systems.drop_duplicates().to_file(38.53729004998249)

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[15], line 1
----> 1 coastal_systems.drop_duplicates().to_file(38.53729004998249)

File ~/mambaforge/envs/codebook/lib/python3.10/site-packages/geopandas/geodataframe.py:1203, in GeoDataFrame.to_file(self, filename, driver, schema, index, **kwargs)
   1120 """Write the ``GeoDataFrame`` to a file.
   1121 
   1122 By default, an ESRI shapefile is written, but any OGR data source
   (...)
   1199 >>> gdf.to_file('dataframe.shp', mode="a")  # doctest: +SKIP
   1200 """
   1201 from geopandas.io.file import _to_file
-> 1203 _to_file(self, filename, driver, schema, index, **kwargs)

File ~/mambaforge/envs/codebook/lib/python3.10/site-packages/geopandas/io/file.py:535, in _to_file(df, filename, driver, schema, index, mode, crs, engine, **kwargs)
    532     df = df.reset_index(drop=False)
    534 if driver is None:
--> 535     driver = _detect_driver(filename)
    537 if driver == "ESRI Shapefile" and any([len(c) > 10 for c in df.columns.tolist()]):
    538     warnings.warn(
    539         "Column names longer than 10 characters will be truncated when saved to "
    540         "ESRI Shapefile.",
    541         stacklevel=3,
    542     )

File ~/mambaforge/envs/codebook/lib/python3.10/site-packages/geopandas/io/file.py:444, in _detect_driver(path)
    442     pass
    443 try:
--> 444     return _EXTENSION_TO_DRIVER[Path(path).suffix.lower()]
    445 except KeyError:
    446     # Assume it is a shapefile folder for now. In the future,
    447     # will likely raise an exception when the expected
    448     # folder writing behavior is more clearly defined.
    449     return "ESRI Shapefile"

File ~/mambaforge/envs/codebook/lib/python3.10/pathlib.py:960, in Path.__new__(cls, *args, **kwargs)
    958 if cls is Path:
    959     cls = WindowsPath if os.name == 'nt' else PosixPath
--> 960 self = cls._from_parts(args)
    961 if not self._flavour.is_supported:
    962     raise NotImplementedError("cannot instantiate %r on your system"
    963                               % (cls.__name__,))

File ~/mambaforge/envs/codebook/lib/python3.10/pathlib.py:594, in PurePath._from_parts(cls, args)
    589 @classmethod
    590 def _from_parts(cls, args):
    591     # We need to call _parse_args on the instance, so as to get the
    592     # right flavour.
    593     self = object.__new__(cls)
--> 594     drv, root, parts = self._parse_args(args)
    595     self._drv = drv
    596     self._root = root

File ~/mambaforge/envs/codebook/lib/python3.10/pathlib.py:578, in PurePath._parse_args(cls, args)
    576     parts += a._parts
    577 else:
--> 578     a = os.fspath(a)
    579     if isinstance(a, str):
    580         # Force-cast str subclasses to str (issue #21127)
    581         parts.append(str(a))

TypeError: expected str, bytes or os.PathLike object, not float

coastal_systems2

	name	lon	lat	zoom	west	south	east	north	geometry
0	Lagoa de Albofeira	-9.181333	38.510379	13	-9.276924	38.483560	-9.085865	-9.085865	POINT (-9.18133 38.51038)
1	Amazon	-50.158433	0.318995	9	-51.913147	-0.780005	-48.405762	1.417092	POINT (-50.15843 0.31899)
2	Columbia river	-124.039699	46.247702	11	-124.477844	46.057509	-123.600998	46.437384	POINT (-124.03970 46.24770)
3	Dune de Pilat	-1.217360	44.586341	13	-1.326942	44.537388	-1.107731	44.635193	POINT (-1.21736 44.58634)
4	Ebro Delta	0.869526	40.699556	12	0.650253	40.595446	1.088676	40.803675	POINT (0.86953 40.69956)
5	Ganges	90.989565	22.312137	8	87.484131	20.267350	94.498901	24.332082	POINT (90.98956 22.31214)
6	Hangzou bay	121.754241	30.488979	9	120.000916	29.537619	123.508301	31.431007	POINT (121.75424 30.48898)
7	Jiangsu coast (north)	119.829333	34.471656	12	119.610214	34.358459	120.048637	34.584888	POINT (119.82933 34.47166)
8	Jiangsu coast (south)	121.020279	32.638839	11	120.581818	32.407212	121.458664	32.869784	POINT (121.02028 32.63884)
9	Mekong river	106.762698	10.179805	11	106.324310	9.909333	107.201157	10.450000	POINT (106.76270 10.17980)
10	Missisipi	-89.245992	29.151004	11	-89.684143	28.910813	-88.807297	29.390551	POINT (-89.24599 29.15100)
11	Pearl (Zhujiang River estuary),	113.636915	22.755011	12	113.417702	22.628272	113.856125	22.881553	POINT (113.63692 22.75501)
12	Redfish pass	-82.200681	26.553486	12	-82.419777	26.430613	-81.981354	26.676300	POINT (-82.20068 26.55349)
13	Sao Fransisco	-48.713194	-26.312026	11	-49.151459	-26.557822	-48.274612	-26.065419	POINT (-48.71319 -26.31203)
14	Senegal delta	-16.667135	13.766087	11	-17.105713	13.499143	-16.228867	14.032679	POINT (-16.66713 13.76609)
15	Severn estuary	-2.680262	51.595222	11	-3.118744	51.424474	-2.241898	51.765715	POINT (-2.68026 51.59522)
16	St Michel	-1.510683	48.633302	13	-1.620312	48.587850	-1.401100	48.678607	POINT (-1.51068 48.63330)
17	Wadden sea	5.195870	53.300644	10	4.319000	52.970973	6.072693	53.627539	POINT (5.19587 53.30064)
18	Wax Lake	-91.448385	29.498476	13	-91.557999	29.438701	-91.338787	29.558228	POINT (-91.44839 29.49848)
19	Willapa bay	-124.071405	46.682020	12	-124.290733	46.587653	-123.852310	46.776082	POINT (-124.07140 46.68202)
20	Yangtze river	121.878402	31.289810	9	120.124512	30.346806	123.631897	32.224419	POINT (121.87840 31.28981)
21	Yellow river	119.173973	37.754205	11	118.735428	37.536955	119.612274	37.971267	POINT (119.17397 37.75420)

coastal_systems2.explore()